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Abstract: This Living Review updates a previous version |Q which its itself 
an update of a review article |Q. Numerical exploration of the properties of 
singularities could, in principle, yield detailed understanding of their nature in 
physically realistic cases. Examples of numerical investigations into the forma- 
tion of naked singularities, critical behavior in collapse, passage through the 
Cauchy horizon, chaos of the Mixmaster singularity, and singularities in spa- 
tially inhomogeneous cosmologies are discussed. 
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1 Introduction 



The singularity theorems |245: 139, 140, 141 1 state that Einstein's equations 



will not evolve generic, regular initial data arbitrarily far into the future or the 
past. An obstruction such as infinite curvature or the termination of geodesies 
will always arise to stop the evolution somewhere. The simplest, physically 
relevant solutions representing for example a homogeneous, isotropic universe 
(Friedmann- Robertson- Walker (FRW)) or a spherically symmetric black hole 
(Schwarzschild) contain space-like infinite curvature singularities. Although, in 
principle, the presence of a singularity could lead to unpredictable measurements 
for a physically realistic observer, this does not happen for these two solutions. 
The surface of last scattering of the cosmic microwave background in the cos- 
mological case and the event horizon in the black hole (BH) case effectively 
hide the singularity from present day, external observers. The extent to which 
this "hidden" singularity is generic and the types of singularities that appear 
in generic spacetimes remain major open questions in general relativity. The 
questions arise quickly since other exact solutions to Einstein's equations have 
singularities which are quite different from those described above. For example, 
the charged BH (Reissner-Nordstrom solution) has a time-like singularity. It 
also contains a Cauchy horizon (CH) marking the boundary of predictability of 
space-like initial data supplied outside the BH. A test observer can pass through 
the CH to another region of the extended spacetime. More general cosmologies 
can exhibit singularity behavior different from that in FRW. The Big Bang in 
FRW is classified as an asymptotically velocity term dominated (AVTD) singu- 



larity 1 94, |165[ since any spatial curvature term in the Hamiltonian constraint 
becomes negligible compared to the square of the expansion rate as the singu- 
larity is approached. However, some anisotropic, homogeneous models exhibit 



Mixmaster dynamics (MD) p3, 187 and are not AVTD — the influence of the 
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spatial scalar curvature can never be neglected. For more rigorous discussions 
of the classification and properties of the types of singularities see |9^, 240 



Once the simplest, exactly solvable models are left behind, understanding 
of the singularity becomes more difficult. There has been significant analytic 



progress 1 244 , 191, 218, ^. However, until recently such methods have yielded ei- 
ther detailed knowledge of unrealistic, simplified (usually by symmetries) space- 
times or powerful, general results that do not contain details. To overcome these 
limitations, one might consider numerical methods to evolve realistic spacetimes 
to the point where the properties of the singularity may be identified. Of course, 
most of the effort in numerica l rela tivity applied to BH collisions has addressed 
the avoidance of singularities |100|. One wishes to keep the computational grid 
in the observable region outside the horizon. Much less computational effort 
has focused on the nature of the singularity itself. Numerical calculations, even 
more than analytic ones, require finite values for all quantities. Ideally then, 
one must describe the singularity by the asymptotic non-singular approach to 
it. A numerical method which can follow the evolution into this asymptotic 
regime will then yield information about the singularity. Since the numerical 
study must begin with a particular set of initial data, the results can never have 
the force of mathematical proof. One may hope, however, that such studies 
will provide an understanding of the "phenomenology" of singularities that will 
eventually guide and motivate rigorous results. Some examples of the interplay 
between analytic and numerical results and methods will be given here. 

In the following, we shall consider examples of numerical study of singulari- 
ties both for asymptotically flat (AF) spacetimes and for cosmological models. 
These examples have been chosen to illustrate primarily numerical studies whose 
focus is the nature of the singularity itself. In the AF context, we shall consider 
two questions: 

The first is whether or not naked singularities exist for realistic matter 
sources. One approach has been to explore highly non-spherical collapse looking 
for spindle or pancake singularities. If the form ation of an event horizon requires 
a limit on the aspect ratio of the matter |239|, such configurations may yield a 
naked singularity. Analytic results suggest that one must go beyond the failure 
to ob serve an apparent horizon to conclude that a naked singularity has formed 
[ 244 1 ■ Another approach is to probe the limits between initial configurations 
which lead to black holes and those which yield no singularity at all (i.e. flat 
spacetime plus radiation) to explore the singularity as the BH mass goes to zero. 
This quest led naturally to the discovery of critical behavior in the collapse of 
a scalar field [0. In the initial study, the critical (Choptuik) solution is a zero 
mass naked singularity (visible from null infinity). It is a counterexample to 
the cosmic censorship conjecture [135]. However, it is a non-generic one since 
fine-tuning of the initial data is required to produce this critical solution. In 
a possibly related study, Christodoulou has shown that for the spherically 
symmetric Einstein-scalar field equations, there always exists a perturbation 
that will convert a solution with a naked singularity (but of a different class 
from Choptuik's) to one with a black hole. Reviews of critical phenomena in 



gravitational collapse can be found in M, 125, 129, 131|. 
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The second question which is now beginning to yield to numerical attack 
involves the stability of the Cauchy horizon in charged or rotating black holes. 
It has been conjectured p45| , [73[ that a real observer, as opposed to a test mass, 
cannot pass through the CH since realistic perturbed spacetimes will convert 
the CH to a strong spacelike singularity [240|. Numerical studies ]5^, p^, |63[ 
show that a weak, null singularity forms first as had been predicted [ pl2| , p02[ [! 

In cosmology, we shall consider both the behavior of the Mixmaster model 
and the issue of whether or not its properties are applicable to generic cosmologi- 
cal singularities. Although numerical evolution of the Mixmaster equations has a 
long history, developments in the past decade were motivated by inconsistencies 
between the known sensitivity to initial conditions and standard measures of the 
chaos usually associated with such behavior [193, 223, 225| , 27, 102 32, 147, 216 1. 
A coordinate invariant characterization of Mixmaster chaos has been formulated 
which, while criticized in its details |194], has essentially resolved the ques- 
tion. In addition, a new extremely fast and accurate algorithm for Mixmaster 



simulations has been developed |37 




We shall discuss results 
14^, p. The key claim 



Belinskii, Khalatnikov, and Lifshitz (BKL) long ago claimed ||17|, |18|, |19|, 
^ that it is possible to formulate the generic cosmological solution to Einstein's 
equations near the singularity as a Mixmaster universe at every spatial point. 
While others have questioned the validity of this claim , numerical evidence 
has been obtained for oscillatory behavior in the approach to the singularity of 
spatially inhomogeneous cosmologies ]25C| , ^ 
from a numerical program to address this issue 
by BKL is that sufficiently close to the singularity, each spatial point evolves 
as a separate universe — most generally of the Mixmaster type. For this to be 
correct, the dynamical influence of spatial derivatives (embodying communica- 
tion between spatial points) must be overwhelmed by the time dependence of 
the local dynamics. In the past few years, numerical simulations of collapsing, 
spatially inhomogeneous cosmolo gical spacetimes have provided strong support 
for the BKL picture 0, || p|]25C|, p|, ||, H). In addition, the Method of 
Consistent Potentials (MCP) [123, |43(| has been developed to explain how the 
BKL asymptotic state arises during collapse. New asymptotic methods have 
been used to prove that open sets exist with BKL's l ocal beha vior (although 
these are AVTD rather than of the Mixmaster type) 163 , 172, |^. Recently, 
van Elst, Uggla, and Wainwright developed a dynamical systems approach to 
G2 cosmologies (i.e. those with 2 spatial symmetries) l242| . 



2 Singularities in AF spacetimes 

While I have divided this topic into three subsections, there is considerable 
overlap. The primary questions can be formulated as the cosmic censorship 



conjecture. The weak cosmic censorship conjecture [209[ requires a singularity 
formed from regular, asymptotically flat initial data to be hidden from an ex- 
ternal observer by an event horizon. An excellent review of the meaning and 



status of weak cosmic censorship has been given by Wald [244[. Counter exam 
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pies have been known for a long time but tend to be dismissed as unrealistic 
in some way. The strong form of the cosmic censorship conjecture |210| forbids 
timelike singularities, even within black holes. 



2.1 Naked singularities and the hoop conjecture 
2.1.1 Overview 

Perhaps, the first numerical approach to study the cosmic censorship conjecture 
consisted of attempts to create naked singularities. Many of these studies were 



motivated by Thome's "hoop conjecture" 239 that collapse will yield a black 



hole only if a mass M is compressed to a region with circumference C < inM 



in all directions. (As is discussed by Wald |244|, one runs into difficulties in 
any attempt to formulate the conjecture precisely. For example, how does one 
define C and M, especially if the initial data are not at least axially symmetric? 
Schoen and Yau defined the size of an arbitrarily shaped mass distribution in 



[228|. A non-rigorous prescription was used in a numerical study by Chiba ]75|.) 
If the hoop conjecture is true, naked singularities may form if collapse can yield 
C > AttM in some direction. The existence of a naked singularity is inferred 
from the absence of an apparent horizon (AH) which can be identified locally. 
Although a definitive identification of a naked singularity requires the event 
horizon (EH) to be proven to be absent, to identify an EH requires knowledge 



of the entire spacetime. While one finds an AH within an EH [166, 167|, it is 



possible to construct a spacetime slicing which has no AH even though an EH 



is present 1 246 1 . Methods to find an EH in a numerically determined spacetime 
have only recently become available and have not been applied to this issue 
[ 179 , 184 . A local prescription, applicable numerically, to identify an "isolated 
horizon" is under development by Ashtckar et al (see for example [Q). 



Figure 1: Heuristic illustration of the hoop conjecture. 



2.1.2 Naked Spindle Singularities? 

In the be st kn own attempt to produce naked singularities, Shapiro and Teukol- 



sky (ST) [230[ co nsid ered collapse of prolate spheroids of collisionless gas. (Naka- 



mura and Sato [ 195 [ had previously studied the collapse of non-rotating de- 
formed stars with an initial large reduction of internal energy and apparently 
found spindle or pancake singularities in extreme cases.) ST solved the general 
relativistic Vlasov equation for the particles along with Einstein's equations for 
the gravitational field. They then searched each spatial slice for trapped sur- 
faces. If no trapped surfaces were found, they concluded that there was no AH 
in that slice. The curvature invariant / — R^upaR^'^'^'^ was also computed. They 
found that an AH (and presumably a BH) formed if C < AttM < 1 everywhere 
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but no AH (and presumably a naked singularity) in the opposite case. In the 
latter case, the evolution (not surprisingly) could not proceed past the moment 



of formation of the singularity. In a subsequent study, ST 231 also showed 
that a small amount of rotation (counter rotating particles with no net angu- 
lar momentum) does not prevent the formation of a naked spindle singularity. 
However, Wald and Iyer p46[ have shown that the Schwarzschild solution has a 
time slicing whose evolution approaches arbitrarily close to the singularity with 
no AH in any slice (but, of course, with an EH in the spacetime). This may 
mean that there is a chance that the increasing prolateness found by ST in effect 
changes the slicing to one with no apparent horizon just at the point required 
by the hoop conjecture. While, on the face of it, this seems unlikely. Tod gives 
an example where a trapped surface does not form on a chosen constant time 
slice — but rather different portions form at different times. He argues that a 
numerical simulation might be forced by the singularity to end before the for- 
mation of the trapped surface is complete. Such a trapped surface would not be 



found by the simulations [241|. In response to such a possibility, Shapiro and 



Tcukolsky considered equilibrium sequences of prolate relativistic star clusters 



[232|. The idea is to counter the possibility that an EH might form after the 
time when the simulation must stop. If an equilibrium configuration is non- 
singular, it cannot contain an EH since singularity theorems say that an EH 
implies a singularity. However, a sequence of non-singular equilibria with rising 
I ever closer to the spindle singularity would lend support to the existence of 
a naked spindle singularity since one can approach the singular state without 
formation of an EH. They constructed this sequence and found that the singular 
end points were very similar to their dynamical spindle singularity. Wald be- 
lieves, however, that it is likely that the ST slicing is such that their singularities 
are not naked — a trapped surface is present but has not yet appeared in their 



time slices [244| 



Another numerical study of the hoop conjecture was made by Chiba et al 
[ [76| . Rather than a dynamical collapse model, they searched for AH's in analytic 
initial data for discs, annuli, and rings. Previous studies of this type were done 



by N akamura et al [|l96| with oblate and prolate spheroids and by Wojtkiewicz 
[ ^5l[ with axisymmetric singular lines and rings. The summary of their results 
is that an AH forms if C < AttM < 1.26. (Analytic results due to Barrabes et al 



|o[ H and Tod |241| give similar quantitative results with different initial data 



classes and (possibly) an alternative definition of C.) 

There is strong analytical evidence against the development of spindle sin- 
gularities. It has been shown by Chrusciel and Moncrief that strong cosmic cen- 
sorship holds in AF electrovac solutions which admit a G2 symmetric Cauchy 
surface [Q . The evolutions of these highly nonlinear equations are in fact non- 
singular. 

2.1.3 Recent results 



Garfinkle and Duncan |114| report preliminary results on the collapse of prolate 



configurations of Brill waves |60| . They use their axisymmetric code to explore 
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the conjecture of Abrahams et al that prolate configurations with no AH 
but large / in the initial slice will evolve to form naked singularities. Garfinkle 
and Duncan find that the configurations become less prolate as they evolve 
suggesting that black holes (rather than naked singularities) will form eventually 
from this type of initial data. Similar results have also been found by Hobill 
and Webster p49|. 

Pelath et al p08| set out to generalize previous results |246, 241 1 that for- 
mation of a singularity in a slice with no AH did not indicate the absence of an 
EH. They looked specifically at trapped surfaces in two models of collapsing null 
dust, including the model considered by Barrabes et al They indeed find a 
natural spacetime slicing in which the singularity forms at the poles before the 
outermost marginally trapped surface (OMTS) (which defines the AH) forms at 
the equator. Nonetheless, they also find that whether or not an OMTS forms in 
a slice closely (or at least more closely than one would expect if there were no 
relevance to the hoop conjecture) follows the predictions of the hoop conjecture. 



Figure 2: This figure is based on Fig. 1 of |208 . The vertical axis is time. The 



blue curve shows the singularity and the red curve the outermost marginally 
trapped surface. Note that the singularity forms at the poles (indicated by 
the blue arrow) before the outermost marginally trapped surface forms at the 
equator (indicated by the red arrow). 



2.1.4 Going further 



Motivated by ST's results [230|, Echeverria [M numerically studied the proper- 



ties of the naked singularity that is known to form in the collapse of an infinite, 
cylindrical dust shell [239|. While the asymptotic state can be found analyti- 
cally, the approach to it must be followed numerically. The analytic asymptotic 
solution can be matched to the numerical one (which cannot be followed all 
the way to the collapse) to show that the singularity is strong (an observer 
experiences infinite stretching parallel to the symmetry axis and squeezing per- 
pendicular to the symmetry axis). A burst of gravitational radiation emitted 
just prior to the formation of the singularity stretches and squeezes in opposite 
directions to the singularity. This result for dust confiicts with rigorously non- 
singular solutions for the electrovac case . One wonders then if dust collapse 
gives any information about singularities of the gravitational field. 

One useful result from dust collapse has been the study of gravitational waves 
which might be associated with the formation of a naked singularity. Such a 



program has been carried out by Harada, Iguchi, and Nakao |16C, 161, 198, 137 



138, 159| 



Nakamura et al (NSN) |197] conjectured that even if naked spindle singu- 
larities could exist, they would either disappear or become black holes. This 
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demise of the naked singularity would be caused by the back reaction of the grav- 
itational waves emitted by it. While NSN proposed a numerical test of their 
conjecture, they believed it to be beyond the current generation of computer 
technology. 



Chiba 1 75 extended previous results to search for AH's in spacetimes 
without axisymmetry but with a discrete symmetry. The discrete symmetry is 
used to identify an analog of a symmetry axis to allow a prescription for an 
analog of the circumference. Given this construction, it is possible to formu- 
late the hoop conjecture in this case and to explore its validity in numerically 
constructed momentarily static spacetimes. Explicit application was made to 
multiple black holes distributed along a ring. It was found that, if the quantity, 
C, defined as the circumference is less than approximately 1.168, a common 
apparent horizon surrounds the multi-black-hole configuration. 

The results of all these searches for naked spindle singularities are contro- 
versial but could be resolved if the presence or absence of the EH could be 
determined. One could demonstrate numerically whether or not Wald's view of 
ST's results is correct by using existing EH finders | 179| , 184| in a relevant sim- 



ulation. Of course, this could only be effective if the simulation covered enough 
of the spacetime to include (part of) the EH. 

2.2 Critical behavior in collapse 
2.2.1 Gravitational collapse simulations 



For a more detailed discussion of critical behavior see [131|. Since Gundlach's 
Living Review article has appeared, the updates in this section will be restricted 
to results I find especially interesting. 

Critical behavior was originally found by Choptuik in a numerical study 
of the collapse of a spherically symmetric massless scalar field. For recent reviews 



see [125, 129 . We note that this is the first completely new phenomenon in 
general relativity to be discovered by numerical simulation. In collapse of a 
scalar field, essentially two things can happen: either a black hole (BH) forms 
or the scalar waves pass through each other and disperse. Choptuik discovered 
that for any 1-parameter set of initial data labeled by p, there is a critical value 
p* such that p > p* yields a BH. He found 

Mbh ^CF{p~P*y (1) 

where Mbh is the mass of the eventual BH. The constant Cp depends on 
the parameter of the initial data that is selected but 7 « .37 is the same for 
all choices. Furthermore, in terms of logarithmic variables p — lnr-|-K, r = 
ln(TQ* — Tg) -|- K (Tq is the proper time of an observer at r = 0, where r is the 
radial coordinate, with Tq* the finite proper time at which the critical evolution 
concludes and k is a constant which scales r) , the waveform X repeats (echoes) 
at intervals of A in T if p is rescaled to /9— A, i.e. X(p — A, r— A) « X(p, r). The 
scaling behavior (|l|) demonstrates that the minimum BH mass (for bosons) is 
zero. The critical solution itself is a counter-example to cosmic censorship (since 



8 



the formation of the zero mass BH causes high curvature regions to become 
visible at r = oo). (See, e.g., the discussion in Hirschmann and Eardley | |145| .) 
The numerical demonstration of this feature of the critical solution was provided 



by Hamade and Stewart |135 . This result caused Hawking to pay off a bet to 
Preskill and Thorne |170t . 



Figure 3: This figure is the final frame of an animation of Type II critical 
behavior in Einstein- Yang-Mills collapse. Note the echoing in the near-critical 
solution. For the entire movie and related references see iTSl. 



Soon after this discovery, scaling and critical phenomena were found in a 
variety of contexts. Abrahams and Evans |^ discovered the same phenomenon 
in axisymmetric gravitational wave collapse with a different value of A and, 
to within numerical error, the same value of 7. (Note that the rescaling of 
r with e'^ ~ 30 required Choptuik to use adaptive mesh refinement (AMR) 
to distinguish subsequent echoes. Abrahams and Evans' smaller A (e^ « 1. 



allowed them to see echoing with their 2-1-1 code without AMR.) Garfinkle [109| 



confirmed Choptuik's results with a completely different algorithm that does 



not require AMR. He used Goldwirth and Piran's |119| method of simulating 
Christodoulou's formulation of the spherically symmetric scalar field in null 
coordinates. This formulation allowed the grid to be automatically rescaled by 
choosing the edge of the grid to be the null ray that just hits the central observer 
at the end of the critical evolution. (Missing points of null rays that cross the 
central observer's world line are replaced by interpolation between those that 
remain.) Hamade and Stewart |135| have also repeated Choptuik's calculation 
using null coordinates and AMR. They are able to achieve greater accuracy and 
find 7 « .374. 

2.2.2 Critical solutions as an eigenvalue problem 

Evans and Coleman |9^ realized that self-similar rather than self-periodic col- 
lapse might be more tractable both numerically (since ODE's rather than PDE's 
are involved) and analytically. They discovered that a collapsing radiation fluid 
had that desirable property. (Note that self-similarity (homothetic motion) is 
incompatible with AF |108| . However, most of the action occurs in the center 
so that a match of the self-similar inner region to an outer A F one sho uld always 



be possible.) In a series of papers, Hirschmann and Eardley |144, 145|1 developed 
a (numerical) self-similar solution to the spherically symmetric complex scalar 
field equations. These are ODE's with too many boundary conditions causing 
a solution to exist only for certain fixed values of A. Numerical solution of this 
eigenvalue problem allows very accurate determination of A. The self-similarity 
also allows accurate calculation of 7 as follows: The critical p = p* solution is 
unstable to a small change in p. At any time t (where t < is increasing toward 
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zero), the amplitude a of the perturbation exhibits power law growth: 



a oc {p-p*)(-ty 



(2) 



where /t > 0. At any fixed t, larger a implies larger Mbh- Equivalently, 
any fixed amplitude a — S will be reached faster for larger eventual Mbh- 
Scaling arguments give the dependence of Mbh on the time at which any fixed 
amplitude is reached: 

Mbh « i-h) (3) 



where 



Thus 



Mbh oc (p - 



(4) 



(5) 



Therefore, one need only identify the growth rate of the unstable mode to obtain 
an accurate value of 7 = 1/k. It is not necessary to undertake the entire dynam- 
ical evolution or probe the space of initial data. (While this procedure allowed 
Hirschmann and Eardley to obtain 7 = 0.387106 for the complex scalar field 
solution, they later found |146| that, in this regime, the complex scalar field has 
3 unstable modes. This means [12J, 131 1 that the Eardley-Hirschmann solution 
is not a critical solution. A perturbation analysis indicates that the critical so- 
lution for the complex scalar field is the discretely self-similar one found for the 



real scalar field |28[.) Koike et al ^7§ obtain 7 = 0.35580192 for the Evans- 
Coleman solution. Although the similarities among the critical exponents 7 in 
the collapse computations suggested a universal value, Maison |183[ used these 
same scaling-perturbation methods to show that 7 depends on the equation of 



state p — kp of the fluid in the Evans-Coleman solution. Gundlach |126| used a 
similar approach to locate Choptuik's critical solution accurately. This is much 
harder due to its discrete self-similarity. He reformulates the model as nonlinear 
hyperbolic boundary value problem with eigenvalue A and finds A = 3.4439. 
As with the self-similar solutions described above, the critical solution is found 
directly without the need to perform a dynamical evolution or explore the space 
of initial data. Hara et al extended the renormalization group approach of [176| 
and applied it to the continuously-self-similar case |l36| . (For an application of 
renormalization group methods to cosmology see [162|.) 



2.2.3 Recent results 



Gundlach |128 completed his eigenvalue analysis of the Choptuik solution to 
find the growth rate of the unstable mode to be 7 = .374 ± .001. He also 
predicted a periodic "wiggle" in the Choptui k m ass scaling relation. This was 
later observed numerically by Hod and Piran | 152 |. Self-similar critical behavior 
has been seen in str ing t heory related axion-dilaton models j9^, 134] and in the 
nonlinear cr-model [ 146 . Garfinkle and Duncan have shown that subcritical 



collapse of a spherically symmetric scalar field yields a scaling relation for the 
maximum curvature observed by the central observer with critical parameters 
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that would be expected on the basis of those found for supercritical collapse 



[1131 



Choptuik et al |79| have generalized the original Einstein-scalar field calcula- 
tions to the Einstein- Yang-Mills (EYM) (for SU{2)) case. Here something new 
was found. Two types of behavior appeared depending on the initial data. In 
Type I, BH formation had a non-zero mass threshold. The critical solution is 
a regular, unstable solution to the EYM equations found previously by Bart- 
nik and McKinnon |0. In Type II collapse, the minimum BH mass is zero 
with the critical solution similar to that of Choptuik (with a different 7 « 0.20, 



A w 0.74). Gundlach has also looked at this case with the same results |127| 



The Type I behavior arises when the collapsing system has a metastable static 



solution in addition to the Choptuik critical one |132|. 

Brady, Chambers and Gongalves [^l], ^ conjectured that addition of a mass 
to the scalar field of the original model would break scale invariance and might 
yield a distinct critical behavior. They found numerically the same Type I and 
II "phases" seen in the EYM case. The Type II solution can be understood 
as perturbations of Choptuik's original model with a small scalar field mass /i. 
Here small means that A/i << 1 where A is the spatial extent of the original 
nonzero field region. (The scalar field is well within the Compton wavelength 
corresponding to /x.) On the other hand, A^ >> 1 yields Type I behavior. The 
minimum mass critical solution is an unstable soliton of the type found by Seidel 



and Suen [229|. The massive scalar field can be treated as collapsing dust to 
yield a criterion for BH formation [ |120| . 

The Choptuik solution has also been found to be the critical solution for 



charged scalar fields ||T3|, |l51[. As p ^ p*, Q/M ^ for the black hole. Q 
obeys a power law scaling. Numerical study of the critical collapse of collision- 
less matter (Einstein- Vlasov equations) has yielded a non-zero minimum BH 



mass 1 217, p01[ . Bizoh and Chmaj |47| have considered the critical collapse of 
skyrmions. 

An astrophysical application of BH critical phenomena has been considered 



by Nimeyer and Jedamzik |199| and Yokoyama [252|. They consider its imph 



cations for primordial BH formation and suggest that it could be important. 
2.2.4 Going further 

The question is then why these critical phenomena should appear in so many col- 
lapsing gravitational systems. The discrete self-similarity of Choptuik's solution 
may be understood as scaling of a limit cycle |136[ . (The self-sim ilari ty of other 



systems may be understood as scaling of a limit point.) Garfinkle [|llO[ originally 
conjectured that the scale invariance of Einstein's equations might provide an 
underlying explanation for the self-similarity and discrete-self-similarity found 
in collapse and proposed a spacetime slicing which might manifestly show this. 



In fact, he later showed (with Meyer) [116| that, while the originally proposed 



slicing failed, a foliation based on maximal slicing did make the scaling manifest. 
These ideas formed the basis of a much more ambitious program by Garfinkle 
and Gundlach to use underlying actual or approximate symmetries to construct 
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coordinate systems for numerical simulations |115[ . 

An interesting "toy model" for general relativity in many contexts has been 
wave maps, also known as nonlinear a models. One of these contexts is crit- 



ical collapse |146 . Recently and independently, Bizoh et al and Liebling 
et al [18C| evolved wave maps from the base space of 3 + 1 Minkowski space 
to the target space . They found critical behavior separating singular and 
nonsingular solutions. For some families of initial data, the critical solution is 
self-similar and is an intermediate attractor. For others, a static solution sep- 
arates the singular and nonsingular solutions. However, the static solution has 
several unstable modes and is therefore not a critical solution in the usual sense. 
Bizoh and Tabor |Q have studied Yang-Mills fields in D -I- 1 dimensions and 
found that generic solutions with sufficiently large initial data blow up in a finite 



time and that the mechanism for blowup depends on D. Husa et al |158| then 
considered the collapse of SU{2) nonhnear sigma models coupled to gravity and 
found a discretely self-similar critical solution for sufficiently large dimensionless 
coupling constant. They also observe that for sufficiently small coupling con- 
stant values, there is a continuously self-similar solution. Interestingly, there is 
an intermediate range of coupling constant where the discrete self-similarity is 



intermittent |238 



Until recently, only Abrahams and Evans [D| had ventured beyond spherical 



symmetry. The first additional departure has been made by Gundlach |130|. 
He considered spherical and non-spherical perturbations of P = perfect fluid 
collapse. Only the original (spherical) growing mode survived. 

Recently, critical phenomena have been explored in 2 + 1 gravity. Pretorius 



and Choptuik |215| numerically evolved circularly symmetric scalar field col- 
lapse in 2 -|- 1 anti-deSitter space. They found a continuously self-similar critical 
solution at the threshold for black hole formation. The BH's which form have 
BTZ exteriors with strong curvature, spacelike singularities in the interior. 
Remarkably, Garfinkle obtained an analytic critical solution by assuming con- 
tinuous self-similarity which agrees quite well with the one obtained numerically 



[tl2| 



2.3 Nature of the singularity in charged or rotating black 
holes 

2.3.1 Overview 

Unlike the simple singularity structure of the Schwarzschild solution, where the 
event horizon encloses a spacelike singularity at r = 0, charged and/or rotating 
BH's have a much richer singularity structure. The extended spacetimes have 
an inner Cauchy horizon (CH) which is the boundary of predictability. To the 



future of the CH lies a timelike (ring) singularity |E45[. Poisson and Israel [212 



213| began an analytic study of the effect of perturbations on the CH. Their goal 
was to check conjectures that the blue-shifted infalling radiation during collapse 
would convert the CH into a true singularity and thus prevent an observer's 
passage into the rest of the extended regions. By including both ingoing and 
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back-scattered outgoing radiation, they find for the Reissner-Nordstrom (RN) 
solution that the mass function (qualitatively Ra/s-yS « M/r^) diverges at the 



CH (mass inflation). However, Ori showed both for RN and Kerr |202, 203 1 that 



the metric perturbations are finite (even though RfivpaR'^"'"^ diverges) so that 
an observer would not be destroyed by tidal forces (the tidal distortion would be 
finite) and could survive passage through the CH. A numerical solution of the 
Einstein-Maxwell-scalar field equations would test these perturbative results. 

For an excellent, brief review of the history of this topic see the introduction 
in |205| . 



2.3.2 Numerical studies 



Gnedin and Gnedin |118t have numerically evolved the spherically symmetric 
Einstein-Maxwell with massless scalar field equations in a 2 -I- 2 formulation. 
The initial conditions place a scalar field on part of the RN event horizon (with 
zero field on the rest). An asymptotically null or spacelike singularity whose 
shape depends on the strength of the initial perturbation replaces the CH. For a 
sufficiently strong perturbation, the singularity is Schwarzschild-like. Although 
they claim to have found that the CH evolved to become a spacelike singularity, 
the diagrams in their paper show at least part of the final singularity to be null 
or asymptotically null in most cases. 

Brady and Smith used the Goldwirth-Piran formulation [119 to study 
the same problem. They assume the spacetime is RN for v < vq. They follow 
the evolution of the CH into a null singularity, demonstrate mass inflation, 
and support (with observed exponential decay of the metric component g) the 
validity of previous analytic results |^12| , 213 , 202 , 203 1 including the "weak" 
nature of the singularity that forms. They find that the observer hits the null 
CH singularity before falling into the curvature singularity at r = 0. Whether 
or not these results are in conflict with Gnedin and Gnedin |118 is unclear Jsof . 
However, it has become clear that Brady and Smith's conclusions are correct. 
The collapse of a scalar field in a charged, spherically symmetric spacetime 
causes an initial RN CH to become a null singularity except at r = where 
it is spacelike. The observer falling into the BH experiences (and potentially 
survives) the weak, null singularity |202, 205 , pl| before the spacelike singularity 
forms. This has been confirmed by Droz ||92| using a plane wave model of the 
interior and by Burko psl using a collapsing scalar field. See also 16^ 



68 



Numerical studies of the interiors of non-Abelian black holes have been car 
ried out by B reitenlohner et al 



7, 58 



104, 105 



and by Gal'tsov et al ||9]j 
(see also |254|1 ) . Although there appear to be conflicts between the two groups, 
the differences can be attributed to misunderstandings of each other's notation 
. The main results include an interesting oscillatory behavior of the metric. 
The current status of the topic of singularities within BH's includes an appar- 
ent conflict between the belief [|9| and numerical evidence |Q that the generic 
singularity is strong, oscillatory, and spacelike and analytic evidence that the 
singularity inside a generic (rotating) BH is weak, oscillatory (but in a different 
way), and null 206 . See the discussion at the end of |206|. 
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Figure 4: Figure 1 from |^6| is a schematic diagram of the singularity structure 
within a spherical charged black hole. 



Various recent perturbative results reinforce the belief that the singularity 
within a "realistic" (i.e. one which results from collapse) black hole is of the 
weak, null type described by Ori [202 , |203[ |. Brady et al |54| performed an anal- 
ysis in the spirit of Bclinskii et al |]19[ to argue that the singularity is of this 
type. They constructed an asymptotic expansion about the CH of a black hole 
formed in gravitational collapse without assuming any symmetry of the per- 
turbed solution. To illustrate their techniques, they also considered a simplified 
"almost" plane symmetric model. Actual plane symmetric models with weak, 



null singularities were constructed by Ori |204|. 

The best numerical evidence for the nature of the singularity in realistic col- 
lapse is Hod and Piran's simulation of the gravitational collapse of a spherically 
symmetric, charged scalar field [154, 153|. Rather than start with (part of) 
a RN spacetime which already has a singularity (as in, e.g., [|5^), they begin 
with a regular spacetime and follow its dynamical evolution. They observe mass 
inflation, the formation of a null sing ularity, and the eventual formation of a 
spacelike singularity. Ori argues |206| that the rotating black hole case is dif- 
ferent and that the spacelike singularity will never form. No numerical studies 
beyond perturbation theory have yet been made for the rotating BH. 

Some insight into the conflict between the cosmological results and those 
from BH interiors may be found by comparing the approach to the singularity 
in Gowdy [121[ spatially inhomogeneous cosmologies (see Sec. |Z|) with T3 
[ ^ and X S'^ [HI I spatial topologies. Early in the collapse, the bound- 
ary conditions associated with the 5^ x topology influence the gravitational 
waveforms. Eventually, however, the local behavior of the two spacetimes be- 
comes qualitatively indistinguishable and characteristic of a (non-oscillatory in 
this case) spacelike singularity. This may be relevant because the BH environ- 
ment imposes effective boundary conditions on the metric just as topology does. 
Unfortunately, no systematic study of the relationship between the cosmological 
and BH interior results yet exists. 



2.3.3 Going further 

Replacing the AF boundary conditions with Schwarzchild-deSitter and RN- 
deSitter BH's was long believed to yield a counterexample to strong cosmic 
censorship. (See [185, 186, 211, ^ and references therein for background and 
extended discussions.) The stability of the CH is related to the decay tails of the 
radiating scalar field. Numerical studies have determined these to be exponen- 
tial [^, ^ rather than power law as in AF spacetimes [^ . The decay tails 
of the radiation are necessary initial data for numerical study of CH stability 
[^ and are crucial to the development of the null singularity. Analytic studies 
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had indicated that the CH is stable in RN-deSitter BH's for a restricted range 
of parameters near extremaUty. However, Brady et al [ p5[ have shown (using 
hnear perturbation theory) that, if one includes the backscattering of outgoing 
modes which originate near the event horizon, the CH is always unstable for all 
ranges of parameters. Thus RN-de Sitter BH's appear not to be a counterexam- 
ple to strong cosmic censorship. Numerical studies are needed to demonstrate 
the existence of a null singularity at the CH in nonlinear evolution. 

Extension of these studies to AF rotating BH's has yielded the surprising 
result that the tails are not necessarily power la w an d differ for different fields. 
Frame dragging effects appear to be responsible |l5C(| . 

As a potentially useful approach to the numerical study of singularities, we 



consider Hiibner's [156, 157 , |l55| numerical scheme to evolve on a conformal 
compactified grid using Friedrich's formalism |103|. He considers the spherically 
symmetric scalar field model in a 2 + 2 formulation. So far, this code has 
been used in this context to locate singularities. It was also used to search for 
Choptuik's scaling |^ and failed to produce agreement with Choptuik's results 
[156|. This was probably due to limitations of the code rather than inherent 
problems with the conformal method. 



3 Singularities in cosmological models 



3.1 Singularities in spatially homogeneous cosmologies 

The generic singularity in spatially homogeneous cosmologies is reasonably well 
understood. The approach to it asymptotically falls into two c lass es. The first, 
called asymptotically velocity term dominated (AVTD) |94, 165 1, refers to a 
cosmology that approaches the Kasner (vacuum, Bianchi I) solution |171 as 
r ^ oo. (Spatially homogeneous universes can be described as a sequence of 
homogeneous spaces labeled by r. Here we shall choose t so that r = oo co- 
incides with the singularity.) An example of such a solution is the vacuum 



Bianchi H model [236| which begins with a fixed set of Kasner-like anisotropic 
expansion rates, and, possibly, makes one change of the rates in a prescribed 
way (Mixmaster-like bounce) and then continues to t = oo as a fixed Kasner 
solution. In contrast are the homogeneous cosmologies which display Mixmaster 
dynamics such as vacuum Bianchi VIII and IX 
and Bianchi I with a magnetic field |178, ^ 
other examples. 




187| , |133| an d Bianchi VIq 
Jantzen |168| ] has discussed 



Mixmaster dynamics describes an approach to the singular- 
ity which is a sequence of Kasner epochs with a prescription, originally due to 
Belinskii, Khalatnikov, and Lifshitz (BKL) |Q, for relating one Kasner epoch 
to the next. Some of the Mixmaster bounces (era changes) display sensitivity 
to initial conditions one usu ally associates with chaos and in fact Mixmaster 
dynamics is chaotic 194 1. Note that we consider an epoch to be a subunit 
of an era. In some of the literature, this usage is reversed. The vacuum Bianchi 
I (Kasner) solution is distinguished from the other Bianchi types in that the 
spatial scalar curvature "^i?, (proportional to) the minisuperspace (MSS) poten- 
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tial 1 187, p27[ , vanishes identically. But arises in other Bianchi types due 
to spatial dependence of the metric in a coordinate basis. Thus an AVTD sin- 
gularity is also characterized as a regime in which terms containing or arising 
from spatial derivatives no longer influence the dynamics. This means that the 
Mixmaster models do not have an AVTD singularity since the influence of the 
spatial derivatives (through the MSS potential) never disappears — there is no 
last bounce. A more general review of numerical cosmology has been given by 
Anninos [§. 

3.2 Numerical Methods 
3.2.1 Symplectic Methods 

Symplectic numerical methods have proven useful in studies of the approach to 



the singularity in cosmological models ||36|. Symplectic ODE and PDE [101, 189 1 
methods comprise a type of operator splitting. An outline of the method (for 
one degree of freedom) follows. Details of the application to the Gowdy model 
(PDE's in one space and one time direction) are given elsewhere po| . 

For a field q{t) and its conjugate momentum p{t) split the Hamiltonian op- 
erator into kinetic and potential energy subhamiltonians. Thus, for an arbitrary 
potential V{q), 

H^\p^ + V{q)^H^{p)+H2{q). (6) 

If the vector X — {p, q) defines the variables at time t, then the time evolution 
is given by 

^ = {H,X}pB^AX (7) 

where { }ps is the Poisson bracket. The usual exponentiation yields an evo- 
lution operator 

for A = Ai + A2 the generator of the time evolution. Higher order accuracy 



may be obtained by a better approximation to the evolution operator [234, 235 [. 
This method is useful when exact solutions for the subhamiltonians are known. 
For the given H, variation of Hi yields the solution 

q = qa+PQAt , p ^ Po, (9) 

while that of H2 yields 

dV 



q = qo , p = PQ- — 
dq 



At . (10) 



Note that H2 is exactly solvable for any potential V no matter how compli- 
cated although the required differenced form of the potential gradient may be 
non-trivial. One evolves from t to t + At using the exact solutions to the sub- 
hamiltonians according to the prescription given by the approximate evolution 
operator (^. Extension to more degrees of freedom and to fields is straightfor- 
ward @, 
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3.2.2 Other Methods 



Symplectic methods can achieve convergence far beyond that of their formal 
accuracy if the full solution is very close to the exact solution from one of the 
subhamiltonians. Examples where this is the case are given in On 
the other hand, because symplectic algorithms are a type of operator splitting, 
suboperators can be subject to instabilities that are suppressed by the full op- 
erator. An example of this may be found in |Q. Other types of PDE solvers 
are more effective for such spacetimes. One currently popular method is itera- 



tive Crank-Nicholson (see [237|) which is, in effect, an implicit method without 
matrix inversion. It was first applied to numerical cosmology by Garfinkle [ [11 1| 
and was recently used in this context to evolve symmetric cosmologies | |44[ |. 

As pointed out in ^ ^ , spiky features in collapsing inhomogeneous 

cosmologies will cause any fixed spatial resolution to become inadequate. Such 
evolutions are therefore candidates for adaptive mesh refinement such as that 



implemented by Hern and Stuart |143, 142| 



3.3 Mixmaster dynamics 
3.3.1 Overview 

Belinskii, Khalatnikov, and Lifshitz (BKL) described the singularity ap- 
proach of vacuum Bianchi IX cosmologies as an infinite sequence of Kasner 



[171 1 epochs whose indices change when the scalar curvature terms in Einstein's 
equations become important. They were able to describe the dynamics approxi- 
mately by a map evolving a discrete set of parameters from one Kasner epoch to 
the next [^, |7^. For example, the Kasner indices for the power law dependence 
of the anisotropic scale factors can be parametrized by a single variable u > 1. 
BKL determined that 

""+^-iK-i)-^ , i<^.„<2 ■ (11) 

The subtraction in the denominator for 1 < m„ < 2 yields the sensitivity to 



initial conditions associated with Mixmaster dynamics (MD). Misner [187| de- 
scribed the same behavior in terms of the model's volume and anisotropic shears. 
A multiple of the scalar curvature acts as an outward moving potential in the 
anisotropy plane. Kasner epochs become straight line trajectories moving out- 
ward along a potential corner while bouncing from one side to the other. A 
change of corner ends a BKL era when u — > {u — Numerical evolution of 

Einstein's equations was used to explore the accuracy of the BKL map as a de- 



scriptor of the dynamics as well as the implications of the map [ 193| , 223, 225, 2^ 



Rendall has studied analytically the validity of the BKL map as an approxima- 



tion to the true trajectories [219| 



Later, the BKL sensitivity to initial conditions was discussed in the language 
of chaos [|ll|, |172| . An extended application of Bernoulli shifts and Farey trees 



was given by Rugh [224| and repeated by Cornish and Levin [M. However, the 



chaotic nature of Mixmaster dynamics was questioned when numerical evolution 
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of the Mixmaster equations yielded zero Lyapunov exponents (LE's) |102, p3, 



147 1 . (The LE measures the divergence of initially nearby trajectories. Only 



an exponential divergence, characteristic of a chaotic system, will yield positive 



exponent.) Other numerical studies yielded positive LE \2U\. This issue was 



resolved when the LE was shown numerically and analytically to depend on the 



choice of time variable |223, Ef|, M. Although MD itself is well-understood, its 



characterization a s ch aotic or not had been quite controversial [148| 



LeBlanc et al [178| have shown (analytically and numerically) that MD can 



arise in Bianchi VIq models with magnetic fields (see also |181|). In essence, the 
magnetic field provides the wall needed to close the potential in a way that yields 
the BKL ma p for u . A similar study of magnetic Bianchi I has been given 
by LeBlanc |l77| . Jantzen ha s dis cussed which vacuum and electromagnetic 
cosmologies could display M D [ll68[ . 

Cornish and Levin (CL) [p6| identified a coordinate invariant way to char- 
acterize MD. Sensitivity to initial conditions can lead to qualitatively distinct 
outcomes from initially nearby trajectories. While the LE measures the expo- 
nential divergence of the trajectories, one can also "color code" the regions of 
initial data space corresponding to particular outcomes. A chaotic system will 
exhibit a fractal pattern in the colors. CL defined the following set of discrete 
outcomes: During a numerical evolution, the BKL parameter u is evaluated 
from the trajectories. The first time u > 7 appears in an approximately Kas- 
ner epoch, the trajectory is examined to see which metric scale factor has the 
largest time derivative. This defines three outcomes and thus three colors for 
initial data space. However, one can easily invent prescriptions other than that 
given by Cornish and Levin |8^ which would yield discrete outcomes. The 
fractal nature of initial data space should be common to all of them. It is not 
clear how the value of the fractal dimension as measured by Cornish and Levin 
would be affected. The CL prescription has been criticized because it requires 
only the early part of a trajectory for implementation JL94[ . Actually, this is 
the greatest strength of the prescription for numerical work. It replaces a single 
representative, infinitely long trajectory by (easier to compute) arbitrarily many 
trajectory fragments. 

Figure 5: The algorithm of is used to generate a Mixmaster trajectory with 
more than 250 bounces. The trajectory is shown in the rescaled anisotropy 
plane with axes /3±/|fi|. The rescaling fixes (asymptotically) the location of the 
bounces. 



To study the CL fractal and ergodic properties of Mixmaster evolution |Q , 
one could take advantage of a new numerical algorithm due to Berger, Garfinkle, 
and Strasser [Q. Symplectic methods are used to allow the known exact solu- 
tion for a single Mixmaster bounce [ ^36| to be used in the ODE solver. Standard 
ODE solvers |214| can take large time steps in the Kasner segments but must 
slow down at the bounce. The new algorithm patches together exactly solved 
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bounces. Tens of orders of magnitude improvement in speed are obtained while 
the accuracy of machine precision is maintained. In |37[ , the new algorithm 
was used to distinguish Bianchi IX and magnetic Bianchi VIq bounces. This 
required an improvement of the BKL map (for parameters other than u) to take 
into account details of the exponential potential. 

So far, most recent effort in MD has focused on diagonal (in the frame of 



the SU{2) 1-forms) Bianchi IX models. Long ago, Ryan |226| showed that off- 



diagonal metric components can contribute additional MSS potentials (e.g. a 



centrifugal wall). This has been further elaborated by Jantzen [ 169 



3.3.2 Recent developments 

The most interesting recent developments in spatially homogeneous Mixmaster 
models have been mathematical. Despite the strong numerical evidence that 
Bianchi IX, etc. models are well-approximated by the BKL map sufficiently 
close to the singularity (see, e.g., p^), there was very little rigorous information 
on the nature of these solutions. Recently, the existence of a strong singularity 
curvatur e blowup) was proved for Bianchi VIII and I X c ollapse by Ringstrom 



[ 221 , |222| and for magnetic Bianchi VIq by Weaver [24J[. A remaining open 
question is how closely an actual Mixmaster evolution is approximated by a 
single BKL sequence [{2191 , |222|| . Since the Berger et al algorithm [^ achieves 
machine level accuracy, it can be used to collect numerical evidence on this 
topic. For example, it has been shown that a given Mixmaster trajectory ceases 
to track the corresponding sequence of integers obtained from the BKL map 
( prl ) at the point where there have been enough era-ending (mixing) bounces to 
lose all the information encoded in finite precision initial data [^ . 

3.3.3 Going further 

There are also numerical studies of Mixmaster dynamics in other theories of 
gravity. For example, Carretero-Gonzalez et al [^ find evidence of chaotic be- 
havior in Bianchi IX-Brans-Dicke solutions while Cotsakis et al [|8^ have shown 
that Bianchi IX models in 4th order gravity theories have stable non-chaotic 
solutions. Barrow and Levin find evidence of chaos in Bianchi IX Einstein- 
Yang-Mills (EYM) cosmologies |lj]. Their analysis may be applicable to the 
corresponding EYM black hole interior solutions. Bianchi I models with string- 
theory-inspired matter fields have been found by Damour and Henneau [89|] to 
have an oscillatory singularity. This is interesting because many other examples 
exist where matter fields and/or higher dimensions suppress such oscillations 
(see e.g. [^). Recently, Coley has considered Bianchi IX brane- world models 
and found them not to be chaotic [Q . 

Finally, we remark on a successful application of numerical Regge calculus 
in 3 -I- 1 dimensions. Gentle and Miller have been able to evolve the Kasner 
solution [117[. 
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3.4 Inhomogeneous Cosmologies 
3.4.1 Overview 

BKL have conjectured that one should expect a generic singularity in spatially 
inhomogeneous cosmologies to be locally of the Mixmaster type (local Mixmas- 
ter dynamics (LMD)) p2| . For a review of homogeneous co smo logies, inho- 
mogeneous cosmologies, and the relation between them, see |l82| . The main 
difficulty with the acceptance of this conjecture has been the controversy over 
whether the required time slicing can be constructed glob ally |]l3 , 122 1 . Montani 
[192 1, Belinskii and Kirillov and Kochnev [175, 174 1 have pointed out that 
if the BKL conjecture is correct, the spatial structure of the singularity could 
become extremely complicated as bounces occur at different locations at differ- 
ent times. BKL seem to imply |22 that LMD should only be expected to occur 
in completely general spacetimes with no spatial symmetries. However, LMD is 
actually possible in any spatially inhomogeneous cosmology with a local MSS 
with a "closed" potential (in the sense of the standard triangular potentials 
of Bianchi VIII and IX). This closure may be provided by spatial curvature, 
matter fields, or rotation. A class of cosmological models which appear to have 
local MD are vacuum universes on x R with a U{1) symmetry [190]. Even 
simpler plane symmetric Gowdy cosmologies [121, |5) have "open" local MSS 
potentials. However, these models are interesting in their own right since they 
have been conjectured to possess an AVTD singularity 1 123 . One way to obtain 



these Gowdy models is to allow spatial dependence in one direction in Bianchi 
I homogeneous cosmologies |2^ . It is well-known that addition of matter terms 
to homogeneous Bianchi I, Bianchi VIq, and other AVTD models can yield Mix- 



master behavior [168, 17S, 177|. Allowing spatial dependence in one direction in 



such models might then yield a spacetime with LMD. Application of this proce- 



dure to magnetic Bianchi VIq models yields magnetic Gowdy models [25C, 247 [ 



Of course, Gow dy cosmologies are not the most general T symmetric vacuum 
spacetimes [121, Restoring the "twists" introduces a centrifugal wall to 

close the MS S. Magne tic Gowdy and general symmetric models appear to 
admit LMD [H^ ^ ~ 



The past few years have seen the development of strong numerical evidence in 
support of the BKL claims [0. The Method of Consistent Potentials (MCP) 



[123 1 has been used to organize the data obtained in simulations of spatially 
inhomogeneous cosmologies [^, |250| , |l[ The main idea is to 

obtain a Kasner-like velocity term dominated (VTD) solution at every spatial 
point by solving Einstein's equations truncated by removing all terms containing 
spatial derivatives. If the spacetime really is AVTD, all the neglected terms will 
be subdominant (exponentially small in variables where the VTD solution is 
linear in the time r) when the VTD solution is substituted back into them. For 
the MCP to successfully predict whether or not the spacetime is AVTD, the 
dynamics of the full solution must be dominated at (almost) every spatial point 
by the VTD solution behavior. Surprisingly, MCP predictions have proved 
valid in numerical simulations of cosmological spacetimes with one [p4| and 
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two ^ Q spatial symmetries. In the case of U{1) symmetric models, a 
comparison between the observed behavior pl| and that in a vacuum, diagonal 
Bianchi IX model written in terms of variables provides strong support for 



LMD 1 45 since the phenomenology of the inhomogeneous cosmologies can be 
reproduced by this rewriting of the standard Bianchi IX MD. 

Polarized plane symmetric cosmologies have been evolved numerically using 
standard techniques by Anninos, Centrella, and Matzner ||^, |6|. The long-term 
project involving Berger, Garfinkle, and Moncrief and their students to study 
the generic cosmological singularity numerically has been reviewed in detail 
elsewhere [^6[ ^ ^ but will be discussed briefly here. 

3.4.2 Gowdy Cosmologies and Their Generalizations 

The Gowdy model on x R is described by gravitational wave amplitudes 
P{9,t) and Q{9,t) which propagate in a spatially inhomogeneous background 
universe described by \{9,t). (We note that the physical behavior of a Gowdy 
spacetime can be computed from the effect of the metric evolution on a test 
cylinder We impose < < 27r and periodic boundary conditions. The 

time variable r measures the area in the symmetry plane with t = oo a curvature 
singularity. 

Einstein's equations split into two groups. The first is nonlinearly coupled 
wave equations for dynamical variables P and Q (where ,a = d/da) obtained 



from the variation of 188 



2tt 



e TTn 



H = ^ J de [n% 



+ ^Jd0 (P,2 + e'^Q,l)] =H,+H, (12) 



where np and ttq are canonically conjugate to P and Q respectively. This 
Hamiltonian has the form required by the symplectic scheme. If the model is, 
in fact, AVTD, the approximation in the symplectic numerical scheme should 
become more accurate as the singularity is approached. The second group of 
Einstein equations contains the Hamiltonian and 0-momentum constraints re- 
spectively. These can be expressed as first order equations for A in terms of P 
and Q. This break into dynamical and constraint equations removes two of the 
most problematical areas of numerical relativity from this model — the initial 
value problem and numerical preservation of the constraints. 

For the special case of the polarized Gowdy model {Q = 0), P satisfies a 
linear wave equation whose exact solution is well-known [ psf . For this case, 
it has been proven that the singularity is AVTD |165[. This has also been 
conjectured to be true for generic Gowdy models U2 



AVTD behavior is defined in |165| as follows: Solve the Gowdy wave equa- 



tions neglecting all terms containing spatial derivatives. This yields the VTD 
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solution [|o| . If the approach to the singularity is AVTD, the full solution comes 
arbitrarily close to a VTD solution at each spatial point as r — > oo. As t cx), 
the VTD solution becomes 



Pi9,T) v{9)t 



(13) 



where t; > 0. Substitution in the wave equations shows that this behavior is 
consistent with asymptotic ex pone ntial decay of all terms containing spatial 
derivatives only if < u < 1 123 . We have shown that, except at isolated 
spatial points, the nonlinear terms in the wave equation for P drive v into this 
range |]35| , |36|| . The exceptional points occur when coefhcients of the nonlinear 
terms vanish and are responsible for the growth of spiky features seen in the 
wave forms We conclude that generic Gowdy cosmologies have an 

AVTD singularity except at isola ted spatial points ||35|, This has been 



confirmed by Hern and Stuart |143| and by van Putten |243f| . After the nature of 
the solutions became clear through numerical experiments, it became possible to 
use Fuchsian asymptotic methods to prove that Gowdy solutions with < t; < 1 
and AVTD behavior almost everywhere are generic |173[ . These methods have 
been recently been applied to Gowdy spacetimes with S'^ x 5*^ and topologies 
with similar conclusions |233|. 

One striking property of the Gowdy models are the development of "spiky 
features" at isolated spatial points where the coefficient of a local "potential 
term" vanishes |3^ . Recently, Rendall and Weaver have shown analytically 
how to generate such spikes from a Gowdy solution without spikes |22C]. 

Addition of a magnetic field to the vacuum Gowdy models (plus a topology 
change) which yields the inhomogeneous generalization of magnetic Bianchi VIq 
models provides an additional potential which grows exponentially if < w < 1 . 
Local M ixmaster behavior has recently been observed in these magnetic Gowdy 
models H^, H^. 

Garfinkle has used a vacuum Gowdy model with S"^ x spatial topology 
to test an algorithm for axis regularity [111|. Along the way, he has shown 
that these models are also AVTD with behavior at generic spatial points that 
is eventually identical to that in the case. Comparison of the two models 
illustrates that topology or other global or boundary conditions are important 
early in the simulation but become irrelevant as the singularity is approached. 

Gowdy spacetimes are not the most general symmetric vacuum cosmolo- 
gies. Certain off-diagonal metric components (th e twists which are gea, ges in 
the notation of (12)) have been set to zero [121|. Restoring these terms (see 
|39|| ) yields spacetim es t hat are not AVTD but rather appear to exhibit a 
novel type of LMD 249 1. The LMD in these models is an inhomogeneous 
generalization of non-diagonal Bianchi models with "centrifugal" MSS potential 
walls [227, 16£| in addition to the usual curvature walls. In [Q, remarkable 
quantitative agreement is found between predictions of the MCP and numerical 
simulation of the full Einstein equations. A version of the code with AMR has 
been developed [^. (Asymptotic methods have been used to prove that the 
polarized version of these spacetimes have AVTD solutions [163 .) 
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3.4.3 U{1) Symmetric Cosmologies 

Moncrief has shown ]19C|] that cosmological models on x R with a spatial 



U{1) symmetry can be described by five degrees of freedom {x, z, A, ip^ui} and 
their respective conjugate momenta {PxjPzjPAiP, All variables are functions 
of spatial variables u, v and time, r. Einstein's equations can be obtained by 
variation of 



H = (p (p dudvH 

dudv {\pI + \e^'pl + y + ie^^r^ - \pl + 2pa) 



- e 



+ e^[(e-^^),„x,„-(e-2^) 



= H1+H2 . (14) 

Here and a; are analogous to P and Q while is a conformal factor for 
the metric eab{x, z) in the u-v plane perpendicular to the symmetry direction. 
Symplectic methods are still easily applicable. Note particularly that Hi con- 
tains two copies of the Gowdy Hi plus a free particle term and is thus exactly 
solvable. The potential term H2 is very complicated. However, it still contains 
no momenta so its equations are trivially exactly solvable. However, issues of 



spatial differencing are problematic. (Curre ntly , a scheme due to Norton [ 200 | 



is used. A spectral evaluation of derivatives [llOO[ which has been shown to work 
in Gowdy simulations does not appear to be helpful in the U{1) case.) A 
particular solution to the initial value problem is used since the general solution 
is not available pQl. 



Figure 6: Behavior of the gravitational wave amplitude at a typical spatial point 



in a collapsing U{1) symmetric cosmology. For details see |41 



Current limitations of the U{\) code do not affect simulations for the po- 
larized case since problematic spiky features do not develop. Polarized models 



have r = ti> = 0. Grubisic and Moncrief |124| have conjectured that these polar 



ized models are AVTD. The numerical simulations provide strong support for 
this conjecture pq, B2I. Asymptotic methods have been used to prove that an 



open set of AVTD solutions exist for this case [164 
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3.4.4 Going further 

The MCP indicates that the term containing gradients of uj in ( p^ ) acts as 
a Mixmaster-hke potential to drive the system away from AVTD behavior in 
generic U{1) models [^^. Numerical simulations provide support for this sugges- 
tion ^ . Whether this potential term grows or decays depends on a function 
of the field momenta. This in turn is restricted by the Hamiltonian constraint. 
However, failure to enforce the constraints can cause an erroneous relationship 
among the momenta to yield qualitatively wrong behavior. There is numerical 
evidence that this error tends to suppress Mixmaster-like behavior leading to 
apparent AVTD behavior in extended spatial regions of U{1) symmetric cos- 
mologies |3^. In fact, it has been found that when the Hamiltonian 
constraint is enforced at every time step, the predicted local oscillatory behavior 
of the approach to the singularity is observed. (The momentum constraint is 
not enforced.) (Note that in a numerical study of vacuum Bianchi IX homo- 
geneous cosmologies, Zardecki obtained a spurious enhancement of Mixmaster 



oscillations due to constraint violation [253, 147 . In this case, the constraint 
violation introduced negative energy.) 

Mixmaster simulations with the new algorithm js^ can easily evolve more 
than 250 bounces reaching \U\ « 10^^. This compares to earlier simulations 
yielding 30 or so bounces with \ « 10*. The larger number of bounces quickly 
reveals that it is necessary to enforce the Hamiltonian constraint. An expli citly 
constraint enforcing [/(I) code was developed some years ago by Ove (see [207| 
and references therein) . 

It is well known that a scalar field can suppress Mixmaster oscillations 
in homogeneous cosmologies. BKL argued that the suppression would also oc- 
cur in spatially inhomogeneous models. This was demonstrated numerically for 
magnetic Gowdy and U{1) symmetric spacetimes [|3|. Andersson and Rendall 
proved that completely general cosmological (spatially T^) spacetimes (no sym- 
metries) with sufficiently strong scalar fields have generic AVTD solutions 
Garfinkle |107| has constructed a 3D harmonic code which, so far, has found 
AVTD solutions with a scalar field present. Work on generic vacuum models is 
in progress. 

Cosmological models inspired by string theory contain higher derivative cur- 
vature terms and exotic matter fields. Damour and Henneaux have applied 
the BKL approach to such models and conclude that their approach to the 
singularity exhibits LMD [|8j . 

Finally, there has been a study of the relationship between the "long wave- 
length approximation" and the BKL analyses by Deruelle and Langlois flOOt]. 



4 Discussion 



Numerical investigation of singularities provides an arena for the close coupling 
of analytic and numerical techniques. The references provided here contain 
many examples of analytic results guided by numerical results and numerical 
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studies to demonstrate whether or not analytic conjectures are vahd. 

Even more striking is the convergence of the separate topics of this review. 
While the search for naked singularities in the collapse of highly prolate sys- 
tems has yielded controversial results, a naked singularity was discovered in the 
collapse of spherically symmetric scalar fields. The numerical exploration of 
cosmological singularities has yielded strong evidence that the asymptotic be- 
havior is local — each spatial point evolves toward the singularity as a separate 
universe. This means that conclusions from these studies should be relevant in 
any generic collapse. This area of research then should begin to overlap with 
the studies of black hole interiors (see for example [Q). 
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